function  results = simulation_other_outputs( sim_matrix, model_sol, par)
% Usage: simulation_other_outputs( sim_matrix, model_sol, par)
% 
% Important:  sim_matrix = [ w_vec,   lambda_vec,   K_vec ];
% Agents take the behavioral lambda and the w as given, but uses the decision rules from the rational model. 

w_vec = sim_matrix(:,1); 
lambda_vec = sim_matrix(:,2); 
K_vec = sim_matrix(:,3);
wb_vec = sim_matrix(:,4); 

%% Parameters
N_w = par.N_w;    N_lambda = par.N_lambda; 
w_grid = par.w_grid;   lambda_grid=par.lambda_grid;  
HQ_Cutoff = par.HQ_Cutoff;

w_dim = kron(ones(N_lambda,1), w_grid );
lambda_dim = kron(lambda_grid, ones(N_w, 1) );
fit_a_function = @(values_matrix)  scatteredInterpolant(w_dim,    lambda_dim ,  reshape( values_matrix , N_w*N_lambda, 1 ),   'linear',   'linear')  ;   % note: linear extrapolation at the boundary

% Portfolio choices
xK_fun =  fit_a_function(model_sol.xK_matrix);
yK_fun =  fit_a_function(model_sol.yK_matrix);
psi_fun = fit_a_function(model_sol.psi_matrix);
rK_fun = fit_a_function(model_sol.muR_matrix + par.AH./model_sol.p_matrix - model_sol.lambda_matrix .* model_sol.kappap_matrix );

% Prices
price_fun =  fit_a_function(model_sol.p_matrix);

% Productivity
productivity_fun = fit_a_function(model_sol.productivity_matrix);

% volatility 
vol_fun =  fit_a_function(model_sol.sigmap_matrix);

% volatility risk premia
vol_premia_fun = fit_a_function( model_sol.xK_matrix .* (model_sol.sigmap_matrix + par.sigmaK) );

% Sharpe Ratio
sharpe_ratio_fun = fit_a_function(  model_sol.sharpe_ratio_matrix    );

% jumps
kappa_GDP_fun = fit_a_function( model_sol.kappa_output_matrix );

% risk-free rate 
rd_fun = fit_a_function( model_sol.rd_matrix );
rf_fun = fit_a_function( model_sol.rf_matrix );

% jumps in prices 
kappap_fun = fit_a_function( model_sol.kappap_matrix );

%% Outputs
results.w = w_vec;
results.lambda = lambda_vec;
results.K = K_vec;
results.wb = wb_vec;

results.xK = xK_fun(  w_vec, lambda_vec );
results.yK = yK_fun(  w_vec, lambda_vec );
results.psi = psi_fun(  w_vec, lambda_vec );

results.p = price_fun(  w_vec, lambda_vec );

results.productivity = productivity_fun( w_vec,  lambda_vec) ;
results.GDP = results.productivity .* K_vec;   
results.bank_credit = results.w .* (results.xK-1) .* results.p .* results.K;    
results.bank_equity = results.p .* results.K .* results.w;   

results.vol = vol_fun(w_vec, lambda_vec);
results.vol_premia = vol_premia_fun(w_vec, lambda_vec);

results.sharpe_ratio = sharpe_ratio_fun(w_vec, lambda_vec);

results.rf = rf_fun(w_vec, lambda_vec);
results.rd = rd_fun(w_vec, lambda_vec);

results.rK = rK_fun(w_vec, lambda_vec);

% newly added on Jun 7, 22. To include the potential jumps in capital price in dN shocks. 
results.kappap = kappap_fun(w_vec, lambda_vec);

results.investment_rate = par.investment( results.p ); 

if( isfield(par, 'supress_credit_spread_in_simulation' ) &&  par.supress_credit_spread_in_simulation )  % do not directly solve credit spread. 
credit_spread_fun = fit_a_function( model_sol.kappap_matrix );   
else
% Credit spread basis for defining credit spread. 
credit_spread_fun = fit_a_function( model_sol.credit_spread );
end

results.credit_spread = credit_spread_fun( w_vec, lambda_vec );
results.mean_spread_before_normalization = mean(results.credit_spread);
results.credit_spread = results.credit_spread/results.mean_spread_before_normalization ;    % the same as Krishnamurthy and Muir to normalize the credit spread

% Then generate bank equity return due to transition. 
equity_transition_return_fun = fit_a_function(model_sol.muw_transition_component);
results.equity_transition_returns = equity_transition_return_fun( w_vec, lambda_vec );


% % ****** Note: newly added. Define a credit spread index by averaging the past
% values of the credit spread. 
% year_back = 10;    % looking back by 10 years to calculate the credit spread.
% looking_back_period = year_back/par.dt;    
% results.credit_spread_index =  movmean(  results.credit_spread,  [ looking_back_period, 0 ]  )   ;
% results.credit_spread_HQ_index =  movmean(  results.credit_spread_HQ,  [ looking_back_period, 0 ]  )   ;
% results.credit_spread_index = zeros( N_sim, 1 );
% results.credit_spread_HQ_index = zeros( N_sim, 1 );
% results.bank_loan_index = zeros( N_sim, 1 );
% results.bank_loan_long_term_index = zeros( N_sim, 1 );
% capital_value_vec =  results.K.*results.p;
% weight_current = 0;
% for( i = 1:N_sim )
%     index_back = i: -1 : max(1, (i-looking_back_period));
%     results.credit_spread_index(i) = (1-weight_current)*mean(  results.credit_spread(index_back) ) + weight_current*results.credit_spread(i) ;
%     results.credit_spread_HQ_index(i) = (1-weight_current)* mean(  results.credit_spread_HQ(index_back) ) + weight_current*results.credit_spread_HQ(i)  ;
%     results.bank_loan_index(i) = (1-weight_current)*mean( capital_value_vec(i)./capital_value_vec(index_back) ) + weight_current * 1 ;
%     results.bank_loan_long_term_index(i) = (1-weight_current)* mean( capital_value_vec(i)./ capital_value_vec(max(1, i-looking_back_period) ) ) + weight_current * 1 ;
% end






